clear all
ssc install genqreg, replace

global total_iterations=250
set seed 113929

use "$dir\DataSCLB", clear

egen gp = group(Group_blind)
xi i.gp

egen sc = group(School_blind)

keep EL_EGRA_PCA_Index MT_Program CCT_Program _Igp_* L_BLabil Dmiss_BL sc gp

save "${dir}/genqreg_data.dta", replace

preserve

	gen int iteration = 0
	
	forvalues q = 5(5)95{
		
		local qtile=`q'/100

		genqreg EL_EGRA_PCA_Index MT_Program CCT_Program, quantile(`qtile') proneness( _Igp_* L_BLabil Dmiss_BL)

		gen double b_MT_Program_q`q' = _b[MT_Program]
		gen double b_CCT_Program_q`q' = _b[CCT_Program]
		
	}	
	
	reg EL_EGRA_PCA_Index MT_Program CCT_Program _Igp_* L_BLabil Dmiss_BL, cluster(sc)
	gen double b_MT_Program_ols = _b[MT_Program]
	gen double lb_MT_Program_ols = b_MT_Program_ols - (invttail(e(df_r),0.025)*_se[MT_Program]) 
	gen double ub_MT_Program_ols = b_MT_Program_ols + (invttail(e(df_r),0.025)*_se[MT_Program]) 	
	gen double b_CCT_Program_ols = _b[CCT_Program]
	gen double lb_CCT_Program_ols = b_CCT_Program_ols - (invttail(e(df_r),0.025)*_se[CCT_Program]) 
	gen double ub_CCT_Program_ols = b_CCT_Program_ols + (invttail(e(df_r),0.025)*_se[CCT_Program]) 	
	
	keep iteration b_MT_Program_q* b_CCT_Program_q* *_MT_Program_ols *_CCT_Program_ols
	keep if _n==1
	
	save "${dir}/genqreg_results.dta",  replace

restore

timer clear 1
timer on 1
forvalues i = 1/$total_iterations{
	preserve
		
		quietly{
			bsample, cluster(sc) strata(gp)
		
			gen int iteration=`i'
		    
			forvalues q = 5(5)95{
				
				local qtile=`q'/100

				genqreg EL_EGRA_PCA_Index MT_Program CCT_Program, quantile(`qtile') proneness( _Igp_* L_BLabil Dmiss_BL)

				gen double b_MT_Program_q`q' = _b[MT_Program]
				gen double b_CCT_Program_q`q' = _b[CCT_Program]
				
			}
			
			keep if _n==1
			keep iteration b_MT_Program_q* b_CCT_Program_q*

			append using "${dir}/genqreg_results.dta"
			
			save "${dir}/genqreg_results.dta", replace
		}
		if mod(`i',10)==0{
		    di "Iteration `i'"
		}
		
	restore
}

timer off 1
timer list 1


use "${dir}/genqreg_results.dta", clear

forvalues q = 5(5)95{

	_pctile b_MT_Program_q`q', p(2.5 97.5), if iteration!=0
	gen lb_MT_Program_q`q' = r(r1), after(b_MT_Program_q`q')
	gen ub_MT_Program_q`q' = r(r2), after(lb_MT_Program_q`q')
	
	_pctile b_CCT_Program_q`q', p(2.5 97.5), if iteration!=0
	gen lb_CCT_Program_q`q' = r(r1), after(b_CCT_Program_q`q')
	gen ub_CCT_Program_q`q' = r(r2), after(lb_CCT_Program_q`q')
	
}

keep if iteration==0

reshape long b_MT_Program_q lb_MT_Program_q ub_MT_Program_q b_CCT_Program_q lb_CCT_Program_q ub_CCT_Program_q, i(iteration ) j(quantile)

rename *_Program_q *_Program_QTE
rename *_Program_ols *_Program_OLS


save "${dir}/genqreg_results_out.dta", replace

# delimit ; 
twoway (rarea ub_MT_Program_QTE lb_MT_Program_QTE quantile, col(gs14) fc(gs14))
	   (line b_MT_Program_OLS quantile, lp(dash) lc(black)) 
	   (line b_MT_Program_QTE quantile, lc(black) lp(solid)) 
	   (line lb_MT_Program_OLS quantile, col(gs11) fc(gs11) lp(dash))
	   (line ub_MT_Program_OLS quantile, col(gs11) fc(gs11) lp(dash)), 
	   graphr(c(white))
	   sch(s2mono)
	   xtitle("Quantile")
	   legend(order(2 3 1) label(1 "95% CI") label(2 "OLS") label(3 "QTE")
	   symx(*0.3) keyg(*0.5) row(1) colg(*0.5)) 
	   ylabel(0(0.5)3.5, angle(horizontal));
# delimit cr
graph export "$output/AppendixFigure2_FullProgram.pdf", as(pdf) replace 


# delimit ; 
twoway (rarea ub_CCT_Program_QTE lb_CCT_Program_QTE quantile, col(gs14) fc(gs14))
	   (line b_CCT_Program_OLS quantile, lp(dash) lc(black)) 
	   (line b_CCT_Program_QTE quantile, lc(black) lp(solid)) 
	   (line lb_CCT_Program_OLS quantile, col(gs11) fc(gs11) lp(dash))
	   (line ub_CCT_Program_OLS quantile, col(gs11) fc(gs11) lp(dash)), 
	   graphr(c(white))
	   sch(s2mono)
	   xtitle("Quantile")
	   legend(order(2 3 1) label(1 "95% CI") label(2 "OLS") label(3 "QTE")
	   symx(*0.3) keyg(*0.5) row(1) colg(*0.5)) 
	   ylabel(0(0.5)3.5, angle(horizontal));
# delimit cr
graph export "$output/AppendixFigure2_ReducedProgram.pdf", as(pdf) replace 

erase "$dir/genqreg_data.dta"
erase "$dir/genqreg_results.dta"
erase "$dir/genqreg_results_out.dta"